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Recent advancements in the synthesis of anisotropic macromolecules and nanoparticles have 
spurred an immense interest in theoretical and computational studies of self-assembly. The cor¬ 
nerstone of such studies is the role of shape in self-assembly and in inducing complex order. The 
problem of identifying different types of order that can emerge in such systems can, however, be 
challenging. Here, we revisit the problem of quantifying orientational order in systems of building 
blocks with non-trivial rotational symmetries. We first propose a systematic way of constructing ori¬ 
entational coordinates for such symmetric building blocks. We call the arising tensorial coordinates 
strong orientational coordinates (SOCs) as they fully and exclusively specify the orientation of a 
symmetric object. We then use SOCs to describe and quantify local and global orientational order, 
and spatiotemporal orientational correlations in systems of symmetric building blocks. The SOCs 
and the orientational order parameters developed in this work are not only useful in performing and 
analyzing computer simulations of symmetric molecules or particles, but can also be utilized for the 
efficient storage of rotational information in long trajectories of evolving many-body systems. 


I. INTRODUCTION 

It is typically easier to visually detect ’order‘ in a particular structure or pattern, than to describe it mathematically. 
What is more difficult- if not impossible- is to rigorously define what constitutes order, as the task of distinguishing 
an ordered structure from a disordered structure can be subjective at times. However, certain types of order such 
as periodicity can be rigorously defined and characterized. With recent advancements in the synthesis of anisotropic 
particles [1-13], it is now possible to assemble more complex forms of ordered structures [14-23]. The problem 
of identifying, distinguishing and quantifying different types of order is therefore of immense practical interest to 
materials science. 

The types of global order that can arise in many-body systems can be loosely classified into two distinct categories. 
Translational order is realized by the positions of the building blocks in a system. Crystalline [24] and quasicrys¬ 
talline [25] order in atomic systems are examples of stand-alone translational order. Rotational order is, however, 
realized by the orientations of the constituent building blocks. In systems of anisotropic building blocks, there is 
usually a strong coupling between translational and rotational order as the position of each building block is dictated 
by the respective shapes and orientations of its neighbors. However, stand-alone rotational order is possible and can, 
for instance, arise in systems of building blocks with large aspect ratios [19]. Nematic liquid crystals are the most 
notable examples, observed in a variety of systems [26-28]. 

There are several well-established methods for characterizing and quantifying translational order. The most popular 
method is to measure or compute the diffraction image [24], /5(q), which is related to the density profile, p(r), via a 
simple Fourier transform: 


p(q) = [ pir)e-^^'^d^r 

The relationship that exists between the symmetries of p(r) and /5(q) is used for determining the symmetry group of 
the crystal. Diffraction images are, however, not sufficient for quantifying the extent of translational order. Bond order 
parameters proposed by Steinhardt et al [29] are widely used for that purpose [30-37]. The extent of translational 
order is quantified by comparing the relative local arrangements of neighboring molecules with that of the ideal 
crystal. The shape matching algorithm proposed by Keys et al [38, 39] is a generalization of Steinhardt’s bond order 
parameters, and can be used to identify and quantify different types of global and local translational order. 
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Characterizing and quantifying rotational order is, however, more difficult. Historically, orientational order param¬ 
eters are obtained from distribution functions [40]. Since the orientation of an arbitrary object is fully determined 
by one polar angle in and three Euler angles in K^, any form of global orientational order can be mathematically 
represented by a distribution function in terms of these angles. Let /(A) be such a function with A = 0 for and 
A = (0, (j), tjj) for K^. Then f{A)dA is the infinitesimal probability that an an arbitrary particle assumes an orientation 
specified by A. If /(•) is a smooth function of A, it can be expanded using a complete set of orthogonal functions in 
the A space. For instance for d = 2 we have: 


fi0) = (1) 

n 

which is the familiar Fourier series expansion of /(0). For d = 3, we have [40]: 

/(0,</>,V') = E (2) 

l,m,n 

where I?[„„’s are the Wigner matrices that form an orthogonal basis for smooth functions of A. If a structure is 
rotationally isotropic, the coefficients of Eq. (1) and (2) will all be zero for nonzero integers. Therefore, the presence of 
any nontrivial nonzero terms in Eq. (1) and (2) is a symptom of broken rotational symmetry. Such nonzero scalars can 
thus be considered as orientational order parameters (OOPs) and collectively describe the orientational distribution 
of the system. In the presence of nontrivial orientational symmetries in the ordered structure, some additional terms 
in (1) and (2) might vanish for symmetry reasons. For instance, in a uniaxial nematic liquid crystal with a director 
along the z axis, will only depend on 0, and thus all the terms corresponding to that explicitly 

depend of (j) and '0 will vanish. That will correspond to terms with m, n ^ 0. The number of order parameters 
are similarly reduced for other point symmetries [40, 41]. This approach has been successfully used to derive order 
parameters for unixial [40] and biaxial nematics [42] and the cubatic phase [43]. It has also been used for obtaining 
suitable expansions for energetic interactions between symmetric molecules and particles [41, 44]. 

The order parameters obtained from the Wigner expansion of a distribution function can be alternatively represented 
as symmetric traceless tensors of different ranks [45, 46]. Therefore each point symmetry group has a collection 
of traceless tensors as its order parameters. The smallest-rank traceless tensor for that particular group is usually 
described as the order parameter for that symmetry group. A list of such tensors for major three-dimensional rotation 
groups are given in Table. 1 of Ref. [46]. 

This classical approach has, however, its own shortcomings. First of all, the mathematical framework for deriving 
such order parameters has been primarily developed for theoretical studies of liquid crystals- such as the ones given 
in [46]- and not for computational studies. Most importantly, those scalar order parameters only show the extent of 
orientational order and provide little evident information about its geometry. For instance, if the director of a nematic 
liquid crystal is not along the z axis- something that almost always happens in a molecular simulation- Eq. (2) will 
be of little practical utility and applying the symmetry constraints to it will not be trivial. For instance, the terms 
with nonzero m or n will no longer vanish under such circumstances. Although certain procedures are in place for 
extracting geometric information for certain phases such as the uniaxial and biaxial nematics, no general framework 
exists for performing this task for an arbitrary type of orientational order. 

Another problem lies in the inherent degeneracy of polar (and Euler) angles in describing the orientations of 
symmetric objects, i.e., objects with nontrivial rotation groups. We schematically depict this in an example in Fig. 1 
where the orientation of a square in two dimensions is unchanged after a 90° rotation while the value of the polar 
angle describing its orientation changes from 45° to 135°. In the conventional approach, these symmetry restrictions 
are imposed on /(A), henceforth limiting the number of non-vanishing terms in Eqs. (1) and (2). A more useful 
alternative, which is the focus of this work, is to replace A with orientational coordinates that are invariant under 
such symmetry operations. This way, any distribution function in terms of such coordinates will automatically satisfy 
the symmetry of the underlying building block(s). 

The main purpose of this work is to propose a systematic procedure for constructing such non-degenerate orien¬ 
tational coordinates for objects with arbitrary rotational symmetries. The proposed coordinates are tensors with 
ranks depending on the rotational symmetry group of the building block. Such coordinates are bijective, i.e., each 
coordinate is associated with one and only one distinct orientation. They are therefore strong descriptors of order, and 
are henceforth called strong orientational coordinates (SOCs). Orientational distribution functions are then expressed 
in terms of such SOCs and not the degenerate polar (or Euler) angles. Order parameters for orientationally ordered 
structures are accordingly derived as the ensemble averages of certain moments of SOCs. The extent and geometry 
of order is then determined from solving a generalized nonlinear optimization problem. 

Besides characterizing and quantifying rotational order, SOCs can be used for measuring spatial and/or tempo¬ 
ral orientational correlations in computer simulations, something that cannot be easily achieved with non-bijective 
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FIG. 1: Degeneracy of polar angles for describing the orientation of a square in two dimensions. 

coordinates (Section IVB). Therefore, SOCs can complement the existing procedures [23] used for defining spatial 
correlation functions for arbitrarily shaped objects, and thanks to their bijectivity, they can provide more accurate 
information about spatial orientational correlations in liquids, glasses and crystals of high-symmetry building blocks. 
Similarly, such coordinates can be the basis of dehning interaction potentials between anisotropic symmetric objects. 
They can also be used for time-averaging the rotational behavior of a collection of particles and/or the entire system 
(Section IV A). This can be used for efficient storage of rotational information for trajectories of evolving many-body 
systems. 

This paper is organized as follows. We introduced the utilized notations and conventions in Section II. We define 
the notion of strong orientational coordinates in Section III A. In Section IIIB, we use the machinery of strong 
orientational coordinates and group theory to derive rotational coordinates for symmetric objects. We then derive 
the corresponding SOCs for all two- and three-dimensional rotation groups. In Section IIIC, we use distribution 
functions of SOCs to quantify different types of rotational order. We then use this framework to define orientational 
order parameters for several liquid crystalline phases alongside a few numerical examples. In Section IV, we discuss 
further potential applications of SOCs other than the quantification of global orientational order. Finally, Section V 
is reserved for concluding remarks. 

We would like to close this introduction by noting that our approach is purely geometrical, i.e., we are not concerned 
about the physical realizability of the described structures. Instead, our aim is to develop computational tools for 
quantifying order in such structures if they ever emerge in simulations. Finally, we will use the terms ’building block’, 
’particle’ and ’object’ interchangeably, all referring to individual symmetric particles in an evolving many-body system. 


II. NOTATIONS AND DEFINITIONS 

The orientational coordinates derived in this paper are all contra-variant tensors with ranks depending on the 
symmetry of the object. In this context, we denote all scalars, i.e., rank-0 tensors, by Greek letters, e.g., a,/3, etc, all 
vectors, i.e., rank-1 tensors, by small italics, like v,w, etc, and all tensors of rank two or higher by capitalized italics- 
e.g. S', T. Also, we will use the conventional Einstein notation for tensors whenever necessary. 

We denote the r-adic power of a vector u G as u’’, i.e., • • • u®'’. Similarly, the r-adic product 

of r distinct vectors wi, U 2 , • • • ,Vr is denoted by uiU 2 • • • For every set of vectors V = {vi,V 2 , - ‘ ‘ its rank-r 

homogeneous tensor, is defined as: 

n 

jr,(v) := (3) 

P=1 

For two tensors S and T of equal rank, a generalized inner product S © T is defined as the full contraction between 
them: 

Here u is the complex conjugate of u. Note that 0 reduces to the Euclidean inner product for r = 1, i.e., when S and 
T are both vectors. Associated with this inner product, a generalized Frobenius norm of a rank-r tensor is defined as: 

\\S\\^ = {S(DSf^ ( 5 ) 

The Frobenius norm reduces to the familiar Euclidean norm for r = 1 and the matrix Frobenius norm for r = 2. 


III. TENSOR ORDER PARAMETERS- GENERAL CONSTRUCTION 

We derive proper orientational order parameters for arrangements of symmetric particles as follows. In Section III A, 
we develop the machinery of homogeneous tensors as tools of describing the orientation of an arbitrary set of vectors 
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in and define the notion of strong orientational coordinates of a collection of vectors. We then map the orientation 
of a symmetric object to a set of equivalent vectors in Section IIIB and use the SOCs derived in Section III A as 
rotational coordinates of symmetric objects. Finally in Section IIIC we use distribution functions of such SOCs to 
identify and quantify order in orientationally ordered arrangements of symmetric objects. 


A. Homogeneous Tensors 

Let V = {vi,V 2 ,- ■ ■ ,Vp} be a finite set of unit vectors in and let W = {Qv : v G V} with Q S 0{d) an 
orthogonal transformation. A function i^(V) is called a strong orientational coordinate of V if ^(QV) = -^(V) 
implies QV = V. Such a function must be invariant under the permutation of the elements of V, a property satisfied 
by homogeneous tensor forms defined in Eq. (3). With the following chain of theorems, we establish that for any V, 
there exists an even number 2q < p and an odd number 2r — 1 < p so that V = QV if and only if J^q{V) = J^q{QV) 
and Jif 2 r-i{V) = Jlf 2 r-i{QV). 

Lemma 1. Let oi, 02 , • • • , a„, 61 , 62 , ■ • • , S K, Then { 01 , 02 , ■ ■ ■ , a„} = {bi, & 2 , ■ • • , bn} (up to multiplicities) if and 
only if o-i = Er=i b} for I <k<n. 

Proof. The proof is given in Appendix A. □ 

Theorem 2. Let V = {vi,V 2 , ■ ■ ■ ,Vp} be an arrangement of p unit vectors in and let W = {Qv : v G V} where 
Q G 0{d) is an orthogonal transformation. Then V = W if and only if JifufY) = for every k < p. 

Proof. ioT 1 < k < p implies that [Jifk{V) — -J^kiW)] © S' = 0 for any S. Set S = v’f and obtain: 

uf © [jfkY) - M{w)] = [[vjvY - [vjwY] = E [4 - 

i=i 

where pij = vfvj and Y = vfwj. Applying Lemma 1 yields {■Cii,''' = {Tug'' id'ip}- However, for every 

1 < i < p, there exists a fij = 1 and this is only possible if there is a vector wj € VF so that Wj = Vi. □ 

This means that the p—tuple {.^i, • • • , .^p) uniquely specifies the orientation of V. However, we are interested 

in decreasing the number of necessary coordinates, ideally to one. The following two theorems refine the scope of our 
search to two coordinates. 

Theorem 3. Let V and W as defined in Theorem 2 and let JifnY) for some integer r. Then 

^r+ 2 (lF). 

Proof. Let := [^(F) - and observe that: 

S Q lJ^n+2(V) - J^n+2(W)] = (F) - (IF) |! ^ > 0 ( 6 ) 

since © S = © [M({V) — J^(IF)] = v’f © [J^(F) — =^(IF)]. Use ( 6 ) and the Cauchy Schwartz 

inequality to conclude that J^+ 2 (F) 7 ^ .^+ 2 (IF) since ||S'|j^ > 0 . □ 

Corollary 4. Let V and W be as defined in Theorem (2) then there exists integers q,r with 2q < p and 2r — I < p 
so that F = IF if and only if J^ 2 q{V) = J^ 2 q{W) and J^r-i(F) = 

Proof. Let Q G 0{d) so that QV 7 ^ F. According to Theorem 2, there exists s < p so that =^(F) 7 ^ JtfsiQV). Denote 
the smallest such integer with sq and define E = {Q; sq is even} and O = {Q-, sq is odd}. Taking 9=5 maxgg^; sq 
and r = i(l + maxggo sq) completes the proof. In the case of either E or O being empty, the proof is completed by 
taking an arbitrary q and r, respectively. □ 

As will be explained below, this is the smallest number of coordinates that can be obtained for an arbitrary set. 
However, further refinement is possible for the subclasses of sets defined below. 

Definition 1. A set F is even if for every v G V, —v G V and odd if for every v G V, —v ^ F. 

Intuitively, one expects the SOC of an even set to be an even-ranked and the SOC of an odd set to be an odd-ranked 
homogeneous form. This is proven in the following theorems. 
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Corollary 5. Let V and W be as defined in Theorem 2. If V and W are even, then there exists an integer 2q < p 
so that V = QV if and only if J^ 2 q{V) = J^° 2 q{QV). 

Proof. Apply Corollary (4) and note that J^ 2 r+i{y) = 0 for every r. □ 

Lemma 6. Let xi, • • • , j/i, • • • , S K with the property that for every distinct i,j < n, if Xi ^ 0, then Xi+Xj 0. 
Then {xi, ■ ■ ■ ,x„} = {yi, ■ ■ ■ ,?/„} (up to multiplicities) if and only ifJ27=i — 127=1 yi^~^ every k < n. 

Proof. As explained in Proposition 4.1 of Ref. [47], the coefficients of p{z) := 11"=!~ Y\a=i{^ ~ Vi) 

can be written as rational functions of Schur staircase functions that only depend on sums of odd powers of Xi ’s (and 
yfs). □ 

Theorem 7. Let V and W- as defined in Theorem 2- he odd sets, then their exists q < p so that = 

J^q-i{W) if and only if V = W. 

Proof. Let J^ 2 q-i{V) = J^ 2 q-i{W) for every q<p and observe that Q[J^ 2 q-i{V) - Jf 2 q-iiW)] = 127=i(lh7~^ ~ 
= 0 for every i < p. Due to the oddness of P, however, pij’s and ^ij's satisfy the conditions specified in Lemma 6 . 
Therefore, {^ij} = {pij} (up to multiplicities) and V = W. The proof is completed by applying Theorem 3. □ 

Remark 1. Note that a set V that is neither even nor odd, can be partitioned into two nonempty even and odd sets 
V = VgU Vq. The SOC for V will therefore be a pair of an even-ranked homogenous form (for 14) and an odd-ranked 
homogeneous form (for 14). An even-ranked form of V cannot exclusively specify the orientation of V as it will 
be invariant under an inversion that will map 14 to —14- Similarly, an odd-ranked homogeneous form will also be 
insufficient as it will be invariant under the orthogonal transformations that alter 14 while keeping 14 unchanged. 

So far, we have established upper bounds on the tensorial rank of strong orientational coordinates of even and odd 
sets. These upper bounds are generic in the sense that they only depend on the cardinality of the underlying set and 
not its structure. In general, tensor SOCs of smaller ranks might be possible for sets with certain symmetries. As will 
become evident in Section IIIB, however, the bounds proposed here can still be tight for certain types of symmetries. 


B. Symmetric Objects 

The SOCs developed in Section III A are permutation-invariant. They can thus be used for describing the orien¬ 
tations of a symmetric object if a bijection can be established between the orientation of the object and a set of 
equivalent vectors. We will achieve this by generating orbits of suitable vectors under the action of the corresponding 
rotation group. It is necessary to emphasize that the construct outlined in Section III A can be applied to any point 
symmetry group, and our consideration of rotation groups is due to the fact that the orientation of a physical object 
(chiral or achiral) can only change as a result of a rotation. In order to establish the connection, it is necessary to 
re-introduce some standard concepts in group theory. 

A rigid body C is symmetric if there is a non-identity Q € SO{d) that The rotation group of a 

rigid body 1% is given by: 


■■= {Q e SO{d) ■.Q£i = £i} (7) 

For every vector v G its orbit is given by: 

:= {Qu:Qg^5?} (8) 

In a Lie group, an orbit can be an uncountably inhnite set. 

Theorem 8. has the following properties: 

1. partitions into equivalency classes i.e. ffy = ffqy for every u G ^ and Q G 

2. is either even or odd. 

Proof. Here we outline the proof for (2) as (1) directly follows from group properties. If for every u G — u ^ ffy, 
then is odd an (2) follows. Suppose there is some u G with —u G ffy and let w G ffy, then there exist 
Q+, Q-, Q G with u = Q+v, —u = Q-V and w = Qv. We thus have: 

- w = -Qv = -QQ'fu = QQf{-u) = QQ+Q-u 

but QQ’)_Q- G and hence —w G Gv and (2) follows. □ 
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FIG. 2: Cn and D„ symmetries in and R®. (a-b) Ce symmetry in (a) R^ and (b) R®. (c) Dq in R^. The vectors shown in 
green are irreducible symmetric descriptors in (a) and (c) but not in (b). 


All vectors in an orbit are equivalent in the sense that they can be mapped onto one another via rotations that leave 
the orientation of ^ unchanged. Each orbit can therefore be uniquely described by an SOC derived in Section III A 
as it is either an even or an odd set. A single orbit is not necessarily sufficient for describing the orientation of We 
might therefore need to map the orientation of a rigid body to a collection of distinct orbits We call 

such a collection a symmetric descriptor of An irreducible symmetric descriptor of is defined as a symmetric 
descriptor that none of its subsets of constituent orbits is sufficient for describing the orientation of Since each 
orbit in an irreducible symmetric descriptor, is uniquely specified by a strong coordinate derived in III A, the 
orientation of can be uniquely specified by the A^-tuple • • • , -^n)- 

There is no unique way of constructing an irreducible symmetric descriptor for a symmetric object as the size and 
the structure of a given orbit will depend on the vector that generates it. In general, one would prefer orbits with 
fewer vectors as their associated SOCs will be of smaller ranks, and will therefore be easier to compute and store on 
a computer. The cardinality of an orbit is at most equal to the order of the rotation group of However, if 
the initial generating vector is chosen so that it is invariant under certain group operation, a smaller orbit will be 
obtained. Such vectors are equivalent to Wyckoff positions in a space group [48]. In general, the orientation of a 
rigid object in can be uniquely specified by at most d — 1 linearly independent vectors. Note that fewer vectors 
might be needed if the rotation group of is a Lie group. A symmetric descriptor should therefore have the following 
properties: 

• It should have sufficient number of linearly independent vectors. 

• It should not be invariant under any rotation that changes the orientation of 

Here we derive irreducible symmetric descriptors and the associated SOCs for all two- and three-dimensional rotation 
groups. The results are summarized in Table. I. 


1. Trivial Rotation Group 

For a non-symmetric object, = {/} and ffy = {u} for every v S An irreducible symmetric descriptor is 
therefore given by with Vi^s being linearly independent. A strong coordinate is therefore given by a (d — 1)- 

tuple = ivi,V 2 , ■ ■ ■ For a non-symmetric object in two and three dimensions, this will 

correspond, as expected, to one and two linearly independent vectors, respectively. 


2. Cn and Dn 

Cn is the only non-trivial rotation group in and corresponds to the symmetry of a regular n-gon (Fig. 2a). In 
three dimensions, C„ corresponds to the symmetry of a pyramid with a regular n-gonal base (Fig. 2b). Z?„, however, 
corresponds to the symmetry of a prism or a bipyramid with a regular n-gonal basis (Fig. 2c). A characteristic 
orbit of both Cn and Dn- denoted by ffn- is generated by an arbitrary unit vector from within the plane that is 
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perpendicular to the n-fold rotation axis. Such an orbit consists of n vectors. In is an irreducible symmetric 

descriptor of Cn- In however, it is only an irreducible symmetric descriptor of Dn (n > 3). For a Cn object, 
however, is invariant under a 180° rotation around one of its constituent vectors, while such a transformation will 
change the orientation of a Cn object. An irreducible symmetric descriptor of Cn is therefore a union of the ’planar' 
orbit, ffni and a single vector that is parallel to the axis of rotation, e.g., the red vector in Fig. 2c. ffn is also not 
an irreducible symmetric descriptor of an object with D 2 symmetry as it only contains two collinear vectors. In that 
case, an irreducible symmetric descriptor can be constructed as the union of {iz} and where 2 ; is the rotation 
axis. 

In Appendix B, we show that J^n is an SOC for the planar orbit Note that for even n, J^fn, and for odd n, 
J^n-i are guarranteed to be SOCs of i^„. Therefore, the upper bound given in Corollary 5 is tight for even ns. 

Remark 2. C 3 is a subgroup of T, the rotation group of a regular tetrahedron, which is a triangular pyramid. D 4 is 
a subgroup of O, the rotation group of an octahedron, which is a square bipyramid, and a cube, which is a square 
prism. These groups have different SOCs and will be discussed separately. 


3. Coo and Doo 

Coo is the symmetry of a cone with a circular base, while corresponds to the symmetry of a cylinder. For both 
symmetries, two types of orbits are possible. An orbit generated by a vector along the rotation axis will be finite and 
will have one and two elements for Coo and Coo, respectively. All other orbits, however, will be uncountably infinite 
sets comprising of one circle for Coo and one or two parallel circles for Doo- It can be noted that a finite orbit 
is an irreducible symmetric descriptor for both Coo and Doo- The corresponding SOC will therefore be = z and 
= zz for Coo and Coo, respectively. 


Tetrahedral Symmetry T 

T corresponds to the rotation group of a regular tetrahedron (Fig. 3) and has twelve elements. A general orbit 
of T will therefore have twelve elements as well. The high symmetry vectors that connect the centroid of a regular 
tetrahedron to its vertices, depicted in green in Fig. 3, can, however, generate an orbit = {ap}p^i that only has 
four elements, is an irreducible symmetric descriptor of a regular tetrahedron. Therefore, its SOC will also be the 
SOC of a regular tetrahedron. Note that is an odd set, and, according to Theorem 7, will be its 

candidate SOCs. However, note that: 






■A ■ 4 ■■ 


(9) 

( 10 ) 


We prove that M 3 is the SOC for Gj- as follows. Let G-r = and Gt' = {bq}q^i be the corresponding orbits 

for the two tetrahedra T and T' and observe that: 


W'M^^Gt) — M3{Gt')\\^p 



p,q=l 


( 11 ) 


with ^pq = Upbq. Contracting (9) and (10) with 6* and yields: 

Ep=i^P, = 0 9=1,---,4 
eU^p, = I 9 = 1,•••,4 


( 12 ) 


As shown in Appendix C, (11) can only be zero if are the permutations of (1,—|,—|,—|) for every q i.e. if 





FIG. 3: A regular tetrahedron with the rotation group T. The smallest orbit has four elements that are depicted in green. 


5 . Octahedral Symmetry O 


O corresponds to rotational symmetry of a cube and an octahedron (Fig. 4) and has twenty-four elements. Orbits 
containing as few as six elements can, however, be constructed by choosing the high-symmetry vectors connecting 
the center of a regular octahedron to its vertices. These vectors are depicted in green in Fig. 4. Such an orbit is an 
irreducible symmetric descriptor of a cube or an octahedron and has the form 0 c = {'^x,±y,±z\ with x,y and z 
being mutually orthogonal. 

Being an even set. Corollary 5 implies that and are the candidate SOCs for 0c- We, however, have: 

= 2{x^xJ +fy^ +z^zJ) = 25^^ (13) 


The only remaining candidates are thus 00/^ and 00^,- We prove that is an SOC of 0c as follows. Let 0c = 
{±Op}p ^2 0c' = {±^g}g=i be the corresponding orbits for the two cubes C and C with and 

note that: 


||jr4(^c)-^4(^C')llF 


16 


3 - E 

P,9=l 


(14) 


with = apbq. For every q we have X]p=i Cpq = 1- As shown in Appendix D, (14) is zero only if ^pg’s are a 
permutation of (1,0,0) for each q i.e. only if 0 c = 0c'- 


6. Icosahedral Symmetry I 


I corresponds to the rotational symmetry of a regular icosahedron and a regular dodecahedron, and has 60 ele¬ 
ments. Orbits with as few as twelve elements can, however, be constructed from the vectors connecting the center 
of an icosahedron to its twelve vertices (Fig. 5). The corresponding orbit denoted by 0i = {dszai}^^^ is an irre¬ 
ducible symmetric descriptor of an icosahedron (dodecahedron). A prototypical set of such unit a^’s are given by 
{a, ±/3,0), (0, a, ±/3), (±^5,0, a) with a = l/-\/l + (j)^, P = (j )/and </> = (-s/S -I- l)/2. The choice of a^’s are 
arbitrary in the sense that both - 1 - 0 ^ and —Oi are valid choices. Also note that of Qj = ±l/-\/5 for i ^ j. 

As an even set, the SOC of 0i will be amongst 002 , ■ ■ ■ , =^2 according to Corollary 5. However one can show that: 


6 



= 2^aX = 4(5*^- 

p^l 

(15) 


6 , 

= 2 V ^ 

p^l 

(16) 
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FIG. 4: A regular octahedron with the rotation group O. The smallest orbit has six elements that are depicted in green. 


With an approach similar to what was used for the tetrahedral and octahedral symmetries, we prove that is the 
SOC for 01 . Let 0i = {±ap}p^i and 0ii = {±.bq}^^i be two such orbits. We have: 


W0i) - Mh{0r)\\l 


128 E 

p,q=l 


156 



(17) 


with ^pq = [apbq]'^. Contracting and ^^(^ 7 ) with bqVq and bqVqbqbq yields: 


6 

~ 2 

p=i 



(18) 

(19) 


According to Appendix E, (17) can only be zero if ^pg’s are a permutation of (1, |) for each q i.e. if 0i 


01 ,. 


C. Ordered Arrangements and Order Parameters 

As explained in Section 1, orientationally-ordered arrangements of symmetric objects can be mathematically de¬ 
scribed using the distribution functions that have been historically expressed in terms of non-bijective orientational 
coordinates such as polar (or Euler) angles. Here, we express such functions in terms of the SOCs derived in Sec¬ 
tion IIIB, and we quantify the extent of orientational order by computing the ensemble averages of certain moments 
of such SOCs. 

An arrangement of (symmetric) objects is called rotationally isotropic if each particle can take all permissible 
orientations with equal probability. This can be characterized by a uniform distribution function. However in an 
orientationally-ordered arrangement of (symmetric) objects, orientational symmetry is broken and each particle tends 
to preferentially take certain orientations more frequently. Such a preference can be characterized by a non-uniform 
distribution function in terms of an SOC. We therefore define a structure- or a phase- as follows. 

Definition 2. A structure or phase of an object is characterized by a probability density function 
where H stands for all the geometric features needed for the macroscopic characterization 

of the structure. 

It is necessary to emphasize that the notion of a phase in Definition 2 entails only the global orientational charac¬ 
teristics of an arrangement, and should not be confused with the thermodynamic notion of a phase that can entail 
both translational and orientational order. 
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FIG. 5: A regular icosahedron with the rotation group I. The smallest orbit has twelve elements that are depicted in green. 

Example 1. Let be a symmetric object with a normalized irreducible symmetric descriptor containing a single 
equivalence class {N = 1) and let M’m be a strong orientational coordinate constructed from jVsi. The density function 
— ’^m) dchues an arrangement of the object ^ where all objects have the same orientation 

with and fl = M’m- 

We now outline the procedure that can be used for deriving orientational order parameters (OOPs) for a phase. 
Let := be an SOC for M and let be a phase. Also consider a tensorial 

function of with the property that 511(0) = satishes the condition that it is a strong descriptor of O, 

i.e., DJli = 9112 if and only if Oi = O 2 . In other words, 911(0) must be invariant under the transformations that keep 
O unchanged. This is to assure that all distinct geometric instances of a phase are distinguishable by 911. However, 
this condition can be relaxed if one is only interested in certain structural features of a phase. 

Now let , ^ 2 ) • ■ ■ ) be an arrangement of (symmetric) particles and define the experimental order estimator 
as: 


_ 1 "■ 

M = (20) 

n ^' 

i=l 

The more perfect the ordering of .^ 1 , .^27 ■ • • 7 -^n is, the closer will M be to 911. Therefore, one would expect ||M— 911||i? 
to be smaller in more perfect arrangements. Henceforth, one can formulate the problem of identifying the underlying 
geometric features of a phase as: 


. ||M-911n||^ 

™ 11911,-911a||^ 

with 911/ = (M)isotropic- One can define scalar order parameter of the phase as: 


( 21 ) 


Kn = 1 - 


||M-911o.|| 

li®l/-®lo*| 


( 22 ) 


where Q* is the minimizer in (21). Note that «: = 0 for a completely isotropic system while k = 1 when matching is 
perfect. For a certain subclass of phases where ||911n|| is constant, one can simplify (21) to: 


max M © 911n 


(23) 
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TABLE I: Strong Orientational Coordinates of Two- and Three-dimensional Rotation Groups. The ‘optimal SOC’ corresponds 
to the smallest-rank strong rotational coordinate. The ‘upper-bound SOC’, however, corresponds to the upper-bound rank 
predicted in Section III A. The vectors given in the column W correspond to a typical orbit. Also note that a = l/^yl + (/>^,/3 = 


^ with (j) = 

(^+l)/2. 




d 


Prototypes 

W 

Optimal SOC 

Upper bound SOC 

2 


regular n-gon 

Vp (cos sin ^) 

P \ 77. ’ 71. / 

w 

II 

Mil, n even 

M 2 n-\, n odd 

3 


{vp}p=i U {z} 

regular n-gonal pyramid (cos 

2(0, 0,1) 

(JC,2) 

= E;=i Vp 

{.Mil, z), n even 
{M 2 p-i,z), n odd 

3 

T>n, n > 2 

regular n-gonal pyramid 
regular n-gonal prism Vp (cos , sin , O) 

A 

w 

II 

Mil, n even 

M 2 n-i, n odd 

3 

T»2 

rectangular parallelepiped 

{± 2 } U {±y} 

2(0, 0,1) 

1/(0,1,0) 

{zz,yy) 

{zz,yy) 

3 

Coo 

Cone, Hemisphere 

{2}, 2(0, 0,1) 

z 

z 

3 

Coo 

Cylinder 

{±2}, 2(0, 0,1) 

zz 

zz 

3 

T 

Regular tetrahedron 

“7S’ 



3 

O 

Octahedron 

Cube 

{±ap}p=i 

ai(l,0,0) 

02 ( 0 , 1 , 0 ) 
03 ( 0 , 0 , 1 ) 

— 2 dpCipCipi 

Op Mi, 

3 

I 

Regular Icosahedron 
Regular Dodecahedron 

{Toplp—]^ 

ai, 2 (a,±d, 0 ) 
03 , 4 ( 0 , a, ±3) 
a5,6(±/3,0,a) 

~ 2 y ^ dpCLpCLpQipi 

dpdp 

In general if the conditions outlined above 

are established for a given fOt, one can use (21) or (23) to obtain SJln* 


and (22) to calculate the scalar order parameter. We will show the procedure of solving (21) or (23) in multiple 
examples at the end of this section. Before doing so, however, we give explicit formulae for calculating 971/ for the 
M’s that are sums of r-adic products. 

Definition 3. Let vi,V 2 ,‘'' , Um S M'’* be arbitrary vectors and oi, 02 , • • • , am be given nonnegative integers adding 
up to a. A symmetrized tensor polynomial ... (r’l, r' 2 ;' • • j I'm) is defined as the sum of all possible direct 

products of the form Vi^Vi^ ''' '>^ia where exactly oi of i^-’s are one, 02 of f^’s are two, etc. The number of distinct 
terms in such a polynomial is a ) ~ Table II gives the list of symmetrized tensor polynomials 

for a < 4 obtained from the following proposition. 

Proposition 9. The symmetrized tensor polynomial defined in Definition 3 is given by: 


Sai.a2,-,a„(^^lW2,"- Wm) 


^Oi,a2,'" + ,a„ ('^1) ’ ’ ’ ) Wfc+lj ‘ ‘ ‘ , Vm) Ofc — 0 

J2k=l ,afc_i,afe-l,afc + i,--- ,a^ (ui, • • • ,Vm) Oi > 0 for all i 


with g ... q{vi,V 2 , - ■ ■ ,Vm) = 1 - 

Proof. The first assertion follows from the definition. For the second assertion look at each term of 
.a„('^1 j ^2, ■ • • ,Vm) and group them based On their initial multiplier. □ 


Proposition 10. Let t G d = 2,3 be a unit veetor and k G Z-°, then : 
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1- {t^^^^)^sotropic = 0 . 

2 ( 2 p)!( 2 fc — 2 p)! 2 fc / fQj. fi — 2 


- ft; K—rr< 

3. {t'^’') isotropic = -, X! X! 


fe k—m / 
\r 


.,n,k—m—n) 


§ 2 m. 2 n, 2 fc- 2 m- 2 n(a:> 2 ^) V = 3. 


2t J_ 1 2l^ /' 2k ' 

,^0 n—0 \2m,2n,2k — 2m — 2n/ 


where (x, y) and (x, y, z) are orthonormal bases for and respectively. 


Proof. Observe that (t") = {1/Sd-i) Jgdt^dfl where Sd-i is the surface of the d-sphere. Contributions from two 
hemispheres to the integral cancel out for odd n and (1) follows. For n = 2k and d = 2 we have t = cos 9 x + sin 6* y 
and: 




1 

27r 


k 

-^2p,2fc-2p§2p,2/c-2p(^> y) 

P—0 


where l 2 p, 2 k- 2 p = /q^’^ cos^^ 0 sin^^ '^^9d6. The odd terms vanish because / 2 p, 2 <j+i = l 2 p+i, 2 q = l 2 p+i, 2 q+i = 0 and 
(2) follows from Eq. (F4). For d = 3, t{9,(j)) = sin0cost/) x + s'mdsmcj) y + cos0 z and: 


{t^^)j 


1 

47r 


k k—p 

EE l2p,2qJ2k-2p-2q,2p+2q+1^2p,2q,2k-2p-2q{x^ 

P—0 q—0 


(a) (2p)!(2(?)! 4P+'^(p + g)!(2fc — 2p — 25 )!/c! 2fc / ^ 

^ ^ AP+ 1 p\q\{p + q)\ (k - p - qy.{ 2 k + ly. ^2p,2q,2k-2p-2q{X,y, Z) 


(a) follows from (F4) and (F7). Rearranging (24) completes the proof. 


(24) 


□ 


This procedure yields the widely known rotation-invariant isotropic tensors given in Table III, and can be thought 
of an algorithmic way of constructing such isotropic tensors for large values of k. Using this proposition, one can 
thus calculate 911/ for any given tensor that is a sum of n-adic products of unit vectors, including the moments of 
homogeneous tensors defined in this work. One can therefore always subtract 911/ in the definition of M so that ||M||^ 
on its own can be used as a measure of how anisotropic a certain arrangement of particles is. The rest of this section 
is devoted to some examples of how (21) and (23) can be formulated and solved. But before doing so, we outline 
the following useful result that can be used to calculate the expected value of a fc-adic power of a vector that can 
uniformly rotate around a rotation axis. 

Proposition 11. Let z,t G M.^ be unit vectors with t uniformly distributed on the plane perpendicular to z. For a 
vector V = az + j3t, (v^) is given by: 


[k/2] I 


m—0 


fc- 2 ; q 2 ; ( 2 w)!( 2 Z- 2 to)! 

— m)\ 


^^2m,2l-2m,k-2l{^i 2 /) 

with X and y being a pair of orthonormal vectors perpendicular to z. 

Proof. Expand in terms of t and z and use case 2 of Proposition 10 to complete the proof. 


□ 


1. Uniaxial Nematics 

Rodlike molecules or nanoparticles can assemble into a rotationally ordered phase known as the uniaxial nematic 
phase in which the rotation axes of all particles are on average aligned to a common vector called a director [26-28]. 
As explained in Section IIIB3, J^({±Zi}) = 2z‘f is the proper SOC for a rod. As for the uniaxial nematic phase, it 
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TABLE II: List ot §ai,a 2 ,a 3 {x,y, zYs for a < 4 and x,y,z£ 


a 

ffli 

02 

03 

^ai,a2,a^ V’’ ■^) 

0 

0 

0 

0 

1 

1 

1 

0 

0 

X 

1 

0 

1 

0 

y 

1 

0 

0 

1 

Z 

2 

2 

0 

0 

XX 

2 

0 

2 

0 

yy 

2 

0 

0 

2 

zz 

2 

1 

1 

0 

xy + yx 

2 

1 

0 

1 

xz + zx 

2 

0 

1 

1 

yz + zy 

3 

3 

0 

0 

XXX 

3 

0 

3 

0 

yyy 

3 

0 

0 

3 

zzz 

3 

2 

1 

0 

xxy + xyx + yxx 

3 

2 

0 

1 

xxz + xzx + zxx 

3 

1 

2 

0 

xyy + yxy + yyx 

3 

1 

0 

2 

xzz + zxz + zzx 

3 

0 

2 

1 

yyz + yzy + zyy 

3 

0 

1 

2 

yzz + zyz + zzy 

3 

1 

1 

1 

xyz + xzy + yxz + yzx + zxy + zyx 


TABLE III: (l^*^)isotropic for a few values of k. 

^il gjkmn _j_ ^im ^jkln _j_ ^in ^jklm 


Here 5'*!*^* = + 5^^8^’’ + 5'^^5^^ and y®!**"*" = ^ijgkimn _|_ gikgjimn _|_ 


k 

^ /isotropic 

/isotropic 

1 


2 ^ ^ijkl 

1 ^ijkl 

2 1 rpijklmn 

1 rpijklmn 

_ 


is fully specified by it, the director, and henceforth, f2 = {it}. In a perfect nematic phase, all particles will align along 
the same director. This perfectly fits into the class of structures described in Example 1. We therefore have: 

- -< 5 *^ 

2 ^ 2 3 

anj} = it*ii^' - - 6 ^^ 
n 3 

p=i 

Since ||91tQ||^ is constant, we can use the optimization problem (23) which takes the form: 


max vFMu 

subject to vFu = 1 

which can be solved by using Lagrange multipliers: 

.if(it, A) = Mu — X{u^u — 1) 

= 2Mu-2Xu = 0 (M-XI)u = 0 


(25) 
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FIG. 6: Nematic OOP (circles) and packing fraction (squares) as a function of dimensionless pressure for a system of 512 hard 
ellipsoids. 


This implies that u should be an eigenvector of M. The largest eigenvalue of M , Ai maximizes Mu. The scalar 
order parameter is given by: 


^nematic 



I) + (A2 + 5) + (A3 + 1 ) 


(26) 


with Ai > A 2 > A 3 being the eigenvalues of M. This formula penalizes any divergence of A^’s from their ’optimal’ 
values of (2/3,—1/3,—1/3) in a perfect nematics. 

As a numerical example, we conduct isothermal isobaric Monte Carlo simulations of a system of 512 hard ellipsoids 
with an aspect ratio a /6 = 4. These highly elongated building blocks are known to assemble into the nematic liquid 
crystal at moderate pressures [49]. Each MC step consists of 512 trial particle sweeps (translation and rotation 
with equal probability), and one box rescale move, on average. The dimensionless pressure is defined based on the 
shortest axis of an ellipsoid that is chosen as the length scale. Snapshots are stored every 20,000 MC steps and the 
nematic order parameter is computed for each configuration. The calculated nematic OOPs are depicted in Fig. 6 . A 
pronounced increase in OOP is observed at around P* = PV/ksT « 5.9, even though the change in packing fraction 
is not as pronounced. 


2. Cubatic Phase 


At sufficiently large densities, hard cubes can assemble into the cubatic phase, a structure in which all particles 
assume the same orientation, while no long-range translational order exists [50]. As explained in Section IIIB5, 
J^ 4 {{±x, ±y, ±z}) is the SOC of a cube. The cubatic phase can be described by = {±Vp}p^i, with Up’s being three 
orthonormal vectors, and M, DJIq and M can be defined as: 


Mijkl 






3 

p^l 

3 

2 E 

P =1 


^ 3 iV 

/^EE<9^ 

p=l q=l 


UpU^u^ul, 


'^p '-'p '-'p 




p,q p,q p,q 


-1 


Ski ^ s^ksji 


(27) 


Here, Up^g’s correspond to the orthogonal vectors describing the orientation of particle q. Since iHo is constant, the 
following optimization problem can be solved for quantifying cubatic order: 

max M Y.p^i (28) 

subject to = 5pg p, g = 1 , 2 ,3 

The constraints of (28) can, however, be equivalently reformulated as X]p=i UpHp = (Multiply both sides by Vg 
and use the linear independence of Up’s to conclude that VpVg = Spg.) The Lagrangian of the optimization problem is 
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therefore given by: 




3 

E 


3 

v^pV^vl - J 2 
i,i = l 


■ 3 


IP^I 


with ^ij = fj,ji. The derivative of the Lagrangian is given by: 


with: 


dvf, 






= M“"'"vivWn + M 


ski 




<! Q Q 


V^VgVg+M 


Tr7b®^ 


V\vy^ + - IM.sV, 


(a) 


bq - CVq 







Vll M12 M13\ 

bq= 4 


, C = 2 

M2I M22 M23 




VM31 M32 /i33/ 


(29) 


(30) 


Note that (a) follows from the invariance of M under index permutation. 'S7q^ = 0 implies that bq = Cvg. Multiplying 
both sides by Vq and summing over q yields: 


EM^ = C'^VqV^ = C (31) 

q=l 9=1 

Note that C is not symmetric for an arbitrary set of orthogonal u^’s and its symmetry is achieved when u^’s are 
amongst the Karush-Kuhn-Tucker (KKT) solutions of (28). The global maximum of (28) can only be obtained if all 
such KKT solutions are identified. We achieve this by adopting a Newton-Raphson scheme that is described below, 
and by performing a sufficient number of attempts using different initial guesses. First, we define (^st ■= “^{l^st — l^-ts) 
for s ^ t and observe that: 


dC,st 

dv^ 


12 


' - M 




V" V 

^p ^p ^p 


)+4(i,; 


.Stm - bibs 


We can therefore approximate C,st as 


(32) 


3 

Cst ~ Cst,0 T ^ ' (^pCst)o ('^p '^p,o) + ■ ■ ■ (33) 

P=1 

The Newton-Raphson iteration can therefore be carried out by simultaneously solving the three equations given by 
(st = 0 and the six equations ensuring the orthonormality of Vq’s, i.e., J^gVqV^ = . We can, however, decrease the 

number of equations from nine to four by using quaternions, which are widely used in simulations of non-spherical 
particles. Unit quaternions are used to describe rigid-body rotations in three dimensions. Rotating a rigid body ^ 
using a unit quaternion q = (( 71 ,( 72 , 93 , 94 ) maps every vector v° G to v = R{q)v°. The rotation matrix i?(q) is 
given by: 


91 + 92 - 93-94 2(9293-9194) 2(9193-^9294) 

2(9194 + 9293) 9 ?- 9 | + 9 i -94 2(9394-9192) 

2(9294-9193) 2(9192-^9394) 91 - 92 - 93 + 94 . 

One can thus consider zip’s- and Cst’s- as implicit functions q. More specifically we have: 


dCst 

dqu 


dCst 


P^l 


dvp 

dqu 




p=i 


dR{q) o 
dqu 
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FIG. 7: Cubatic OOP (circles) and packing fraction (esquires) vs. dimensionless pressure for a system of 512 hard cubes. 


where v° correspond to a set of orthonormal vectors corresponding to q = 1. In order to solve C = under the 
orthonormality constraint, one can solve the following four equations denoted by f(q) = 0: 


f(q) 


( \ 

Cl 3 

C 23 

\q 1 + «! + + 94 -1/ 


The Newton-Raphson iteration can therefore be carried out using the following formula: 


with: 


Qn+l = q™ - ^(qn) ^f(qn) 


$(q) 


^dCi2ldqi dCi2ldq2 
dCu/dqi 5C13/992 
9 ( 23 /dqi 8 ( 23 /dq 2 

V 2 qi 2q2 


dCi 2 /dq 3 8(12/dqA 
8Ci3/8q3 8Ci3/8q4 
8C,23/dq3 dC,23/8q4 

293 2^4 / 


(34) 


(35) 


Iteration (35) can be carried out for a number of randomly-selected unit quaternions as initial guesses. Once the global 
maximum is attained, Eq (22) can be used to calculate the scalar cubatic order parameter. Neither the accuracy nor 
the convergence rate depends on the selection of as they are only benign parameters of f(q). 

As a numerical example, we conduct isothermal isobaric Monte Carlo simulations of a system of 512 hard cubes in 
the range of pressures at which the cubatic phase is known to exist. Technical specifications of the MC simulations 
are identical to what we discussed earlier for the hard ellipsoid system. The dimensionless pressure is computed using 
the edge length of a cube as the length scale. In order to increase the numerical efficiency of the Newton-Raphson 
iterations, backtracking is also utilized. Fig. 7 depicts packing fraction and cubatic OOP as a function of dimensionless 
pressure. The isotropic-to-cubatic transition occurs at P* = PV/kT « 6.1 and is marked by pronounced jumps in 
both the packing fraction and the cubatic OOP. 


3. Tetratic Phase 

Similar to the cubatic phase, hard squares that are the two-dimensional equivalents of hard cubes can form a 
rotationally ordered phase with C 4 symmetry, known as the tetratic phase, at sufficiently high pressures [51, 52]. As 
explained in Section IIIB2, J^({±a:, ±?/}) is the SOC for the C 4 symmetry. The tetratic phase can be described by 
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FIG. 8: (74 OOP (cirlces) and packing fraction (squares) vs. pressure for a system of 256 hard squares. 


ft = {±Vp}p^i, with vi and V 2 being two orthonormal vectors, and M, 91to and M can be defined as: 


Mnki ^ ^ + 5^^5^^) 

p=i 

2 1 




2 N 




(36) 


p=l q=l 


Here, Up^^’s are the corresponding orthogonal vectors describing the orientation of particle q. In order to identify the 
Wl that best matches M, one can solve the following two-dimensional equivalentt of the optimization problem given 
by Eq. (28): 


max Yll=i 

subject to VpV^ = Spq 


^p ^p 


P: 9 = 1, 2 


(37) 


with the solution given by C 12 = /^12 ~ ^-21 = 0. Numerical estimation of the Karush-Kuhn-Tucker solutions of (37) 
can be performed using a similar Newton-Raphson scheme, based on 9, the polar angle by noting that: 


with: 


C^Cl2 


p=l 


P 


m 


( cos 9 sin 9 
— sin 9 cos 9 


(38) 


(39) 


As a numerical example, we conduct isothermal isobaric Monte Carlo simulations of a system of 256 hard squares 
in the range of pressures at which the tetratic phase emerges The simulation details, including what constitutes an 
MC step, is identical to what was discussed in the hard ellipsoid and the hard cube systems. Dimensionless pressure 
is computed from the edge length of a square as the length scale. Simulations are performed for 10® MC cycle, and 
snapshots are gathered every 5 x 10"^ MC cycles. Fig. 8 depicts the packing fraction and the tetratic OOP as a function 
of pressure. The isotropic-to-transition occurs gradually for P* = PA/kT « 8.0 and is not as pronounced as the 
transitions observed in three-dimensional systems. 


4- Tetrahedral Nematics 

Under external fields, regular tetrahedra are capable of forming a structure in which all particles take the same 
orientation. The possibility and the thermodynamic nature of such a transition has been thoroughly discussed in 
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theoretical works of liquid crystals [46]. However, such a transition has yet to be realized experimentally. The 
arising orientationally ordered structure is typically referred to as the tetrahedral nematics phase, a structure that 
geometrically falls into the general category of structures outlined in Example 1. According to Section IIIB4, 
is the SOC for the tetrahedral rotation group. For the tetrahedral nematic phase, H = {up}p^i is the characteristic 
orbit of the reference tetrahedron with which the individual tetrahedra align. M, Slln and M are thus defined as: 




m 


ijk 

n 




p^l 

4 

p^l 

^ N 3 
1=1 P=1 


(40) 


Here, Vp^q^s are the elements of the high-symmetry orbit for particle q. The following optimization problem is to be 
solved for quantifying the extent of tetrahedral order: 

max M l^p^^ u’pulu^ 

subject to VpVq = ^6pq — ^ p, g = 1,2 ,3,4 

As for the cubatic phase, we can replace the constraints of (41) with X)p=i = (4/3)d*-^ and X)p=i = 0. (See 
Appendix G for details.) The Lagrangian and its derivatives are thus given by: 


4 4 

^ UpU^pUp - Y 

P—1 p—1 

AJ qA£ = bq — C — CUq 

with: 



S 

CO 




^ Til T12 Ail3\ 

bq = 


,c = 

A 2 

,c = 

T21 T22 T23 




{X3) 


\T31 T32 T33) 


with C = Multiplying Vq.5f by on the right and summing over q yields: 


^ (42) 

9=1 

The symmetry of C can be enforced using a method similar to what was explained for the cubatic phase and the 
scalar order parameter can be calculated accordingly. 


5. Tetrahedral Axial Nematics 

We call an arrangement of regular tetrahedra a tetrahedral axial nematic if a specific axis of each particle (modulo 
symmetry operations) aligns with a common director on average. The formation of tetrahedral axial nematics has 
neither been observed in experiments, nor has it been suggested in theoretical studies. The discussion that follows is 
therefore only a geometrical illustration of the ideas presented in this paper. 

Like the uniaxial nematic phase described in Section HIC 1, = {±z} for a tetrahedral axial nematic. Here, we 

consider three distinct plausible types of alignments with z. In the face nematic phase, each tetrahedron has one of 
the vectors that is normal to its faces aligned with z. In the edge nematics, however, one edge of each tetrahedron 
is aligned with z. The third structure is what we call the z nematics, and a vector that connects the centers of two 
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TABLE IV: Projection parameters and ui for face, edge and 2 : uniaxial nematics. 


Phase 

ai 

Pi 

0^2 

P2 

0:3 

Ps 

0:4 

Pi 

V^4 4 

^j=l 

V^4 ^ 2 ,d2 

.^i —1 Pi 


U 

Face Nematics 

Edge Nematics 

z nematics 

1 

3 

V3 

3 

0 

3 

V6 

3 

1 

3 

V6 

3 

V3 

3 

V8 

3 

v/3 

3 

vq 

J2 

1 

3 

0 

V3 

3 

Vs 

3 

1 

V6 

3 

1 

3 

0 

n/3 

3 

Vs 

3 

1 

\/6 

3 

28 

27 

8 

9 

4 

2 

8 

27 

4 

9 

8 

2 

64 

27 

20 

9 

16 

2 

28 

27 

7 

18 

14 

2 


TABLE V: Expected values of and gnir) for different nearest neighbors shells in the hard tetrahedron system 


nth nearest neighbor 


9H{r) 

1 

32 

1 

81 

9 

9 

2464 

77 


2187 

243 

Q 

79072 

2471 


59049 

6561 


non-adjacent edges aligns with 2 . Note that all these phases have inversion symmetry. Therefore, M cannot be an 
odd-ranked moment of Instead we define M as: 








N 4 


- vEE 

9=1 p=i 


vl, „vi 




Note that -I- (4/9)5®-^so M is indeed a function of For each of the three phases 

introduced above, ftRn can be calculated by expressing the elements of the high-symmetry orbit as Vp = 

ttpZ + j5ptp with tp T z, and using Proposition 11 to obtain: 




E A i j k I 
p^l 


a^z z z -\- 


1 


4 ^4 . 

Z 40 


(43) 


p=i 


p=i 


with P and Q given by: 

pijki _ x’' + y^y^y'^y^) + x^'x^y’^y^ + y^y^x’^x^ + x''y^y^x^ + y'^x^x^y’' + x'^y^x^y’^ + y^x^y'^x'' (44) 

Q^jki _ x'^x^z^z^ + y^y^z'^z^ + z^z^x^x^ + z‘z^y^y^ + z^x^x^z^ + z^y^y^z^ 

(45) 


+x'"z^ z^x^ + y'^z^ z^y^ + z'^x^ z^x^ + z'^y^ z'"y'‘ + x^’z^x^z^ + y'^z^y^z'' 


Here, x and y are mutually orthogonal unit vectors that are also perpendicular to z. It can be easily observed that 
P and Q are invariant under orthogonal transformations that keep z unchanged. Therefore, SDlo clearly satisfies the 
properties outlined in Section IIIC. The extent of order can therefore be quantified by solving the following associated 
optimization problem: 


max 




subject to z*z* = 1 


(46) 


where u> = J2p=i[^p l^^p ~ The derivation details can be found in Appendix H. Table IV gives the 

corresponding a^’s, /?p’s and w’s for the face, edge and z nematics. For the z nematic phase, a; < 0, and solving (46) 

reduces to minimizing z"^ z^ z^ z^ under the same constraints. The scalar order parameter can also be defined 

using (22). 


IV. OTHER APPLICATIONS OF STRONG ROTATIONAL COORDINATES 

So far, we have only used SOCs to quantify global orientational order in systems of particles with non-trivial 
rotational symmetries. In this section, other potential applications of SOCs in computational studies of soft condensed 
matter are discussed. 
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FIG. 9: gnir) calculated for a fluid of 4,096 hard tetrahedra at different packing fractions. 



r 


FIG. 10: g{r) and gnir) calculated for a quasicrystal and a fluid of 8,000 hard regular tetrahedra at Pa^/ksT = 64, with a 
the edge length of a tetrahedron. 


A. Time-averaged orientations 

The machinery of strong orientational coordinates can be used to calculate the average orientation of a particle 
in the course of a simulation. In general, this is done to get rid of thermal fluctuations, and can be considered the 
equivalent of determining the inherent structure in an energetic system [53]. More precisely, if Mp{t) is the SOC of 
particle p at time t, one can maximize ^fftp © ^ fj' Mp(t)dt to solve for the best average orientation that matches the 
trajectory of p. This problem also fits into the general class of problems outlined in Example 1. 


B. Spatial Correlation of Local Rotational Order 

There are different ways of quantifying spatial orientational correlations, with different levels of generality [23]. 
The SOCs derived in this work offer a systematic way of quantifying orientational correlations in theoretical and 
computational studies of building blocks with nontrivial symmetries. We explain this through an example, i.e., the 
different phases formed by hard tetrahedra. At low packing fractions, hard tetrahedra form a simple fluid in which the 
orientations of different particles are not correlated. This, however, changes as the packing fraction increases, and a 
networked structure emerges in which all neighboring particles are in face-to-face configurations. At sufficiently high 
densities, this network transforms into a dodecagonal quasicrystal [15]. In both the complex fluid and the quasicrystal, 
the neighboring tetrahedra are in face-to-face contact, a fact quantitatively characterized in earlier studies [23]. The 
SOCs derived in this work can, however, allow us to get further information about the transformation of the simple 
fluid into the complex fluid vis-a-vis the angular distribution of the neighboring face-to-face tetrahedra. It can also 
be used for detecting and quantifying long-range orientational coherence in these networks. 
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It was shown in Section IIIB4 that is an SOC for a regular tetrahedron. Now consider an arrangement of N 

tetrahedra and let N{r) be the number of pairs whose center-to-center distance lies in the interval [r — dr 12^ r + dr/2] 
and define the following correlation function: 


9H{r) 


N{T)\\m\ 


T E 

F <p^q->^\d{p,q)—r\<dr/2 


i k -^✓>7 2 i k 
^3,p ^3,9 


(47) 


We can quantify different types of spatial correlations in the orientations of tetrahedra using (47). We first analyze 
different scenarios for the first nearest neighbor shell. Geometrically, the smallest distance between the centroids of two 
neighboring tetrahedra is achieved when their faces are touching one another. At low densities, the two neighboring 
tetrahedra can freely rotate around the axis perpendicular to this touching face, while at higher densities, rotation 
is restricted and the two faces tend to match perfectly. The value of gnir) can be calculated for these two idealized 
configurations. Let 2 be a unit vector perpendicular to the common plane of the two touching faces. The characteristic 
vectors of each tetrahedron can thus be expressed as: 


Ui = 

z 

Up = 


Vi = 

—z 

Vp = 

3 tp 


with: 

Sp = cos[27r(p — 2)/3]a: + sin[27r(p — 2)/3]j/ 

tp = cos [27r(p — 2)/3 + 0] a: + sin [27r(p — 2)/3 + 0] 2 / 

and X and y two orthonormal vectors perpendicular to z. For a perfect non-rotating face-to-face contact, 9 = 0 and 
= —If • The corresponding value of gnir) will thus be —| « —0.11111. For a freely-rotating face-to-face 
contact, however, 6 ~ [/(0,27r), and which yields a g//(r) value of —| « 0.5555. Investigating 

the magnitude of the first valley of gH{f) will elucidate the type of the face-to-face contact prevalent in that particular 
phase. Fig. 9 shows gH{r) vs. r for three different packing fractions. At </> = 0.24, the first valley of gnix) has a value 
of —0.503 which is very close to the calculated value of — | for a perfect freely rotating face-to-face configuration. At a 
packing fraction of 0.48 however, the first valley of gnix) has a value of —0.15 which is very close to the theoretically 
predicted value of — | for the perfect face-to-face configuration. At intermediate densities (the violet curve), the first 
valley is in between — | and — |. This analysis can be extended to the second and third nearest neighbor shells. For 
perfect face-to-face configurations, the corresponding values of 9 h{x) are given in Table V. The gnix) 

values shown in Fig. 9 are consistent with these theoretical predictions, which confirms the existence of the network 
structure in which all neighbors are in perfect face-to-face contact. 

Differences in spatial orientational correlation can be the basis of differentiating the quasicrystal and the complex 
fluid, both of which have a similar network structure with perfect face-to-face contacts between nearest neighbors [15, 
16]. They also have very similar radial distribution functions, g{x)'s, as depicted in Fig. 10. Note the lack of long- 
range order in both (?(r)’s. To the contrary, these phases have starkly different gnix) functions. For the disordered 
fluid, orientational correlations are virtually nonexistent at large separations. For the quasicrystal, however, such 
correlations never decay and are long-ranged. By using gnix), we are therefore able to detect long-range orientational 
coherence in a networked phase of hard tetrahedra. It is noteworthy that the gnix) values corresponding to the first, 
second and third nearest neighbor shells in the quasicrystal are closer to the theoretically predicted values given in 
Table V than the corresponding peaks in the fluid. Therefore, the face-to-face contacts in the quasicrystal are more 
perfect than the fluid. 


V. CONCLUSION 

In this work, we revisit the problem of quantifying orientational order in arrangements of anisotropic building blocks, 
and we propose a systematic way of constructing symmetry-invariant coordinates for arbitrary building blocks. We 
call the arising tensorial coordinates strong orientational coordinates, and we discuss their potential applications in 
theoretical and computational studies of symmetric building blocks. For instance, we demonstrate that such SOCs 
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can be used in the systematic quantification of spatiotemporal correlations in colloidal systems. Most importantly, 
the orientational distribution functions can be expressed in terms of such SOCs, and the problem of identifying, 
quantifying and describing long-range rotational order can be formulated as a generalized non-linear optimization 
problem. The arising scalar and tensorial order parameters can be efficient ways of storing rotational information on 
a computer, especially for objects with nontrivial symmetries. 
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Appendix A: Proof of Lemma 1 

Form the polynomials p(z) — — ai) and q{z) = YYi^i(,z — bi) and expand them to obtain p{z) = z'^ + 

and q{z) = z^ + where the coefficients are given by Newton’s identities [54]: 

k-l 

(t>k = 

k-l 

where (pk = V'fc = ^i- Note that (pi = tpi if and only if ai = /3i and by induction ak = Pk for every 

1 < k < n. Thus p{z) and q{z) have the same coefficients, and thus the same roots. 


Appendix B: Proof of J^n being an SOC for the planar orbit G-n in Section IIIB2 


We prove this through the following chain of lemmas and theorems. 

Lemma 12. Let vi^V’i^--- ,Vn S be unit vectors so that VpVq = cos [27r(p — gj/nj; then for n > 3, Vi’s span a 
two-dimensional subspace ofM.‘^. 

Proof. Let x,y,z G be unit vectors so that vi = x,V 2 = cos(27r/n)a: -|- sin(27r/n)y and Vk = cos[27r(fc — l)/n]a; -|- 
sin[27r(fc — l)/n]z for some k >3. Use v'^V 2 = cos[27r(fc — 2)/n] to deduce z = y. □ 

Lemma 13. Letp,q be nonnegative integers and (f>k = 0-\-2nk/n. Then Ip{ 6 ) = cos^ sin"^ i/fc is independent 

of 0 for p q < n, and is a function of 0 if p q = n. 

Proof. Note that: 


cos^ (pk 


sin'^ pk 


(0 ^ ~ 



(—1)"* exp [ipki‘ 2 m 


9)] 


(Bl) 

(B2) 


lp,q{0) can therefore be written as: 


lp,q{0) 


i’i 

2P+9 


p q 

EE 


m—O 




(-ire 


i6{2l-\-2m—p—q) 


^ ^ |^^ 27 ri( 2 /+ 2 m—p—g)/n 
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Poj. g277^(2^^-2m-p-g)/r^ ^ 1 have: 


^ |^g 27 ri{ 2 ^+ 2 m—p—g)/n 


_ ^2'n'i{2l-\-2m—p—q) 

^ g 27 r 2 ( 2 Z+ 2 ? 7 i—p—g)/n 


Thus the sum over k survives only if 21 + 2m — p — q is a multiple of n. However \2l + 2m — p — q\ < p + q- Thus if 
p + q < n, the only possibility is zero, which will take away the 0 dependence of lp^g{0). For p + q = n, two 0-dependent 
terms will survive i.e. and the proof follows. □ 


Theorem 14. is an SOC of ffn- 

Proof. Let V = = {vk}k=i- Without loss of generality suppose Vk = cos [9 + 2TTk/n] ex + sin [9 27rfc/n] with 

Cx and Cy being the unit vectors along the x and y directions. The components of are therefore either zero or 
Ip,m-p which are always independent of 9 for m < n according to Lemma 13. For even n, the proof follows from 
Corollary 5. The proof for odd n is completed by noting that Jifn{QV) will have some nonzero components that are 
zero for if QV and V are not in the same plane. □ 


Appendix C: Derivation Details of Tetrahedral SOC 


We need to solve the following optimization problem: 

maximize 

subject to = 1 ^ = 1,2,3,4 (Cl) 

= 0 Z = l,2,3,4 


The form of (Cl) suggests that we can solve it by breaking into four independent problems of the form: 


maximize 
subject to 


ELi Xi = 0 
= 1 


(C2) 


We solve this problem by enumerating all solutions satisfying the Karush-Kuhn-Tucker criteria [55, 56] and show that 
global maximum is attained if a:: is a permutation of (1, ——|). The Lagrangian for (C2) is: 


.if (x, v, jj) 

dxj 


4 4 

+ v^x, + p 


r 4 




4 

3 


3a:^ + V + 2pxj 


KKT conditions requires that Va;.if = 0. Eliminating p one easily obtains: 



i=l 


which implies that: 

-p± ^p'^ +3 


(C3) 

(C4) 


(C5) 


(C6) 


Denote the roots of this equation with and and let n+(n_) be the number of xfs equalling ^_|_(^_). Using 
-I- n_^_ =0 one gets: 




3m 


•^3(16 — m^) 
m =F 4 

-y/3(16 — m^) 


(C7) 

(C8) 


where m = — n_. The global maximum is thus obtained for m = —2,^+ = 1,^_ = — for which Ei=i *? “ §■ 

The global maximum of (Cl) is therefore 
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Appendix D: Derivation Details of Octahedral SOC 


Letting (pq = ^pq, we need to solve the following optimization problem: 

max El,q=l Cpq 

subject to X]p=i Cp 9 = 1 <? = 1,2,3 (Dl) 

Cp 9 > 0 p,g = l,2,3 

Observe that (Dl) can be broken into three independent and yet identical optimization problems of the form: 

E 3 2 

i=i^i 

subject to 2^1 = 1 (D 2 ) 

Xi > 0,i = 1,2,3 


which is solved by identifying x’s that satisfy the Karush-Kuhn-Tucker criteria [55, 56]. The Lagrangian is given by: 


3 3 

2=1 2=1 


—— = 2xj — i/o — A = 0 
dxj 


3 

X'^Xi 

i=l 


(D3) 

(D4) 


which yields Xj = {i'j+X)/2. Let n < 3 be the number of nonzero ccj’s. For each such xj, Vj = 0 due to complementary 
slackness. We thus have A = ^ and Xj = A for nonzero Xj's. Choosing n = 1 maximizes the objective function with 
Si=i Xi = 1. The global maximum for (Dl) is therefore 3 with ^pqS being a permutation of (±1, 0, 0) for each q. 


Appendix E: Derivation Details of Icosahedral SOC 

We need to solve the following optimization problem: 

max EL=i ^p .9 

subject to X]p=i Cp 9 = 2 g = 1, • • • , 6 

y^6 C2 ^ 6 2 • • • 6 

Z^p=i spg 5 y 5^ 

Cpg > 0 = 1, • • • ,6 

which can be broken into six independent and yet identical optimization problems of the form: 

max Y.i=i xt 

subject to XlLi xl = I 

ELi Xi = ‘2 

Xj>0 j = 1, 2, • • • , 6 

In order to solve (E2), we first solve the following optimization problem: 

max E"=i xl 
subject to Yh=i xj = i 
EEi^* = 2 

using Lagrange multipliers. The Lagrangian is given by: 

n 

•^{,x\ , ■ * *, Xfi, X, fi') 'y ] x^ A 

i=l 

which can be differentiated to get: 

—— = 3a:^ — 2 Xxj — u = 0 
dxq ^ 




.2 = 1 


y^^xi -2 


. 2=1 


(El) 


(E2) 


(E3) 


(E4) 
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Summing over j yields M = ^ [^ — 4A] . The roots of the quadratic equation (that we denote by and ^_) are thus 
given by: 


3 




(E5) 


Let n+(n_) be the number of Xj’s being equal to and let m = n+ — ri-. From n+^+ + n_^_ = 2 we have: 


A= ® 
n 


c± = - 

n 


1 — mi. 


' 3n - 10 


1 — (m =F n)4 


3n- 10 


10(?t.2 — m?) 

which suggests that the problem is not feasible for n < 4. The value of the objective function is given by: 

36 16 8m(3n — 10) j 3n — 10 


n+^;+n_C_ = 2 

on 


5n^ 


10(77,2 — m 2 ) 


(E6) 

(E7) 

(E8) 


Since d[m/\/n‘^ — m'^]/dm = j\J (n2 — 7772)2 > 0, the function is maximized for the smallest possible m i.e. m = 2—77 

and the global maximum is given by: 


36 16 4(77 - 2)(3n-10) 3n-10 

577 772 5772 y 10(77 — 1) 


(E9) 


Now we solve (E2) by identifying x’s that satisfy the Karush-Kuhn-Tucker criteria. Its Lagrangian is of the form: 

6 


.if = ^ x^ — ViXi — A 


i=l 

„2 


Ed- 


.7=1 




y^xj -2 




o 2 

gijT — 3Xj — Vj — 2Xxj — fj, 


It can be easily shown that for every Xj = 0, Vj = and is independent of j. For Xj ^ 0, complementary slackness 
implies that Vj = 0 and jdxj = 0 implies fJ, ^ 0- Note that non-zero xj^s are amongst the Karush-Kuhn-Tucker 
solutions of (E3) for 77 = 6 — ttq with uq being the number of Xj’s that are zero. Since /(4) < /(5) < /(6) from 
Eq. (E9), 77o = 0 gives the global maximum with X]i=i ~ ~ value of the objective function 

of (El) is thus ^ which can be only achieved for ^p^’s being permutations of (1, |). 

Appendix F: Some Useful Trigonometric Integrals 

We are interested in integrals of the form: 


r27r 


d-m.n — 


cos’" 0 sin” 0 de 


J„„ = / cos’" 6» sin” 
Jo 


(FI) 

(F2) 


Note that l 2 m-i, 2 n = l 2 m, 2 n-i = 0 since the integrand is odd around 0 = tt. Also: 


l2ra-l.,2n-l = 2 / COS^’”"" 0 sin^""' 0 d0 = 0 

Jo 


since the integrand is odd around 0 = f. I 2 m, 2 n{fn,n > 0) can be calculated by considering the identity: 


72m,0 — 


27r(2777)! 

4’"(m!)2 


(F3) 
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and integration by part: 


cos2™0sin2"6»d6» = 


Cos2m-l0gin2"+l0l2^ 


2 tz “t- 1 


2m — 1 
27Z “t" 1 


^2m — 2,2n+2 — 


+ 1^ / cos2'"”^6'sin2"+2 0d6» 

0 "A 

27r(2m)!(2n)! 

4'"+"'m!n!(m + n)! 


(F4) 


(F3) can be proven as follows: 


J2m,0 — 




27r 12m 

P j 


since exp [2i(p — m)] (i0 = 0 for p 7^ m. Note that J 2 m+i,n = 0 since the integrand is odd around n/2. For 
m = 2p, n = 2q both even we have: 


J'. 


2p,2q 


cos^^ X sin^'* xdx = If 


2p,2q 


(F5) 


We also have: 


J: 


2m ,1 — 


cos^*" X sin xdx = - 


— COS^'"^^ X 


2m + 1 

Similar to what was done for /2m,2n; using integration by part we have: 


2m + 1 


J2m 


2n 


,2n+l — 


2m + 1 


J: 


2n 2n — 2 


2m+2,2n—1 — 


2 2 2^"+^n!(2m)!(m + n)! 

2m + 1 2m + 3 2m + 2n — 1 2m + 2n + 1 m!(2m + 2n + 1)! 


(F6) 


(F7) 


Appendix G: Equivalency of Rigidity constraints for regular tetrahedron 

We need to show that the following set of constraints are equivalent: 


p^l 


(Gla) 


E 


u; = o 


(Gib) 


4 1 

u^pUl = -Spq - -,p,q = 1,2, 3, A 


(G2) 


To prove that (Gl) implies (G2), multiply (Gla) by Uq and observe that '^p^i{u^Uq)up = |ug, which, after some 
rearrangement takes the form: 


E 


2 + ^pq ipp 


Up = 0 


with ^pq = alaq. However since X)p=i WpitJ is a nonsingular matrix, any three of four ttp’s are linearly independent 


and ^qq — Cp9 = I • We also have 


Z)p=i UpUp = 16/3 which yields ^qq = l,^pq = —^{p ^ q) and (G2) follows. 


To prove (Gl) from (G2) take an arbitrary set of vectors satisfying (G2) and show that (Gl) holds; however since 
0 and d*^’s are isotropic tensors, they are invariant under any unitary transformation and the proof follows. 
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Appendix H: Derivation Details of Uniaxial Nematics Order Parameter 


The objective function of the associated optimization problem, can be simplified by observing that M 

is both invariant under index permutation and traceless. More specifically: 




The same thing can be done for the six terms containing z, x and y. With a similar argument one can verify that: 





The objective function can therefore be written as: 
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